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1.  INTRODUCTION 


It  is  axiomatic  that  all  phases  of  future  semiconductor  device  develop¬ 
ment,  such  as  that  indicated  by  the  VLSI  (very  large  scale  integration),  and 
VHSIC  (very  high  speed  integrated  circuits)  programs  depend  on  the  quality  of 
the  semiconductor  material.  Growers  of  these  materials,  e.g.,  silicon,  gallium 
arsenide,  and  indium  phosphide,  etc.,  are  interested  in  efficient,  reliable 
methods  that  yield  defect-free  and  dislocation-free  semiconductors.  Over  many 
years,  a  variety  of  different  techniques  have  been  developed  to  grow  materials, 
some  overcoming  disadvantages  or  shortcomings  of  earlier  methods,  while  others 
are  designed  to  overcome  specific  problems  associated  with  the  growth  of  a  par¬ 
ticular  semiconductor.  Examples  here  include  the  Czochralski  pull  method  which 
was  developed,  among  other  reasons,  to  minimize  temperature  variations  in  the 
melt.  The  Czochralski  method  is  one  of  the  most  common  crystal  growing  techniques 
used  in  silicon  technology  and  single  crystals  three  or  four  inches  in  diameter 
are  common  today.  Another  example  involves  a  modification  of  the  Czochralski 
method  and  incorporates  liquid  encapsulation.  This  is  finding  increased  use  in 
gallium  arsenide  technology,  where  it  tends  to  eliminate  arsenic  in  the  growth 
system.  The  development  and  use  of  these  growth  systems  and  others  (e.g.,  Bridgeman 
zone  refining,  floating  zone,  vapor,  liquid  and  molecular  beam  epitaxy,  chemical 
vapor  deposition,  etc.),  has  been  largely  empirical  due  in  a  major  part  to  the 
absence  of  physical  theories  with  predictive  capabilities.  These  difficulties 
lie  in  two  groups.  The  first  involves  the  incomplete  picture  of  growth  mechanisms. 
The  second,  while  interrelated  to  the  first,  recognizes  the  fact  that  growth  is 
fundamentally  a  three-dimensional  phenomena  and  that  one-  or  two-dimensional 
theoretical  approximations  of  this  phenomena  will  not  necessarily  predict  in  a 
meaningful  way  such  things  as  the  local  temperature  distribution  near  the  solid- 
liquid  interface,  temperature  distributions  within  the  melt  due  to  heating  coils, 
the  distribution  of  unwanted  impurities  arising  from  contamination  by  the 
crucible,  etc . 

The  question  then  becomes:  How  can  one  successfully  undertake  a  program 
that  attempts  to  simulate  crystal  growth  and  thereby  provide  a  basic  physics 
for  this  growth.  An  attempt  as  a  partial  answer  to  this  question  is  the 
subject  of  this  report  which  summarizes  Scientific  Research  Associates  Phase  I  DESAT 


study  for  simulating  Gochralski  crystal  growth.  The  Phase  I  study  relied  strongly 
on  the  attitude  that  fluid  flow  in  the  melt  is  one  of  the  major  influences  on  the 
final  properties  of  the  grown  crystal  [Refs.  1-10].  As  a  consequence,  it  dealt  solely 
with  obtaining  solution  to  the  relevant  hydrodynamic  equations  for  examining  flow 
patterns  in  the  melt. 


2.  PHASE  I  TECHNICAL  OBJECTIVES 

\J 

The  primary  purpose  of  the  Phase  I  study  was  to  establish  the  feasibility 
of  developing  reliable  theories  with  predictive  capability  in  the  area  of  crystal 
growth.  Within  the  guidelines  of  the  DESAT  program  the  proposed  study  was 
limited  to  an  area  where  some  computational  fluid  dynamic  studies  have  already 
appeared  and  where  numerous  experimental  results  have  been  reported;  namely, 
the  area  of  Czochralski  crystal  pulling.  The  intent  of  the  study  was  to 
generalize  all  previously  developed  work,  and  to  incorporate  three-dimensional 
effects,  which  invariably  arise  when  stirring  is  included.  For  Czochralski 
growth  nonaxisymmetric  effects  were  formulated  and  numerical  simulations  were 
performed  for  two  cases:  (1)  a  local  hot  spot,  on  the  crucible  wall,  and 
(2)  angular  misalignment  between  the  crystal  and  crucible  rotational  axes.^ 

A ' 
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3.  NUMERICAL  CONSIDERATIONS 


In  the  Czochralski  method  of  crystal  growth,  the  growing  crystal  is 
pulled  slowly  from  the  center  of  a  rotating  crucible  of  melt  material  and  the 
resultant  crystal  properties  are  strongly  influenced  by  details  of  the  melt 
thermal  and  velocity  fields.  The  melt  flow  represents  a  relatively  complex 
fluid  mechanics  pehnomenon.  As  shown  in  the  schematic  of  Fig.  1,  the  melt  is 
contained  in  a  nominally  axisymmetric  crucible  of  constant  temperature  which 
may  be  rotating  about  its  axis.  The  crystal  which  develops  at  the  melt  upper 
surface  is  nominally  concentric  with  the  crucible  and  both  axes  are  nominally 
aligned.  The  crystal  rotates  about  its  axis  during  the  growth  process.  The 
upper  surface  of  the  melt  is  a  free  surface  whose  shape  is  determined  by  the 
constant  pressure  free  surface  boundary  condition  and  a  meniscus  of  unknown 
shape  occurs  at  the  melt  crystal  boundary.  Heat  transfer  at  the  melt-crystal 
boundary  is  determined  by  heat  conduction  away  from  this  boundary  within  the 
crystal  and  the  heat  transfer  at  the  melt  free  surface  is  determined  by 
radiative  effects  to  the  atmosphere.  Finally,  the  melt  motion  is  driven  both 
by  inertial  and  viscous  forces  associated  with  the  crystal  and  crucible  rotation 
and  by  bouyancy  forces  arising  due  to  temperature  differences  within  the  melt. 
Obviously,  with  all  these  phenomena  affecting  the  flow,  a  prediction  of  the 
melt  velocity  and  temperature  fields  is  a  formidable  task. 

To  date,  several  investigations  have  focused  upon  predictions  of  the 
melt  flow  field.  In  a  relatively  early  work  Kobyashi  and  Arizumi  (Ref.  2) 
considered  axisymmetric  flow  situations  and  solved  the  steady-state  stream 
function-vorticity  form  of  the  Navier-Stokes  equations  in  the  melt  and  a 
heat  conduction  equation  in  the  crystal  to  determine  the  shape  of  the 
crystal-melt  interface.  The  analysis  included  the  effects  of  bouyancy,  but 
assumed  the  melt  free  surface  to  be  flat.  Using  a  numerical  over-relaxation 
technique  to  solve  the  equations,  the  authors  obtained  converged  solutions 
and  examined  the  effects  of  fluid  motion  on  the  interface  shape  for  viscous 
bulk  Czochralski  flows. 
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In  a  more  recent  effort  Langlois  and  Shir  (Ref.  6)  applied  the  unsteady 
Navier-Stokes  equations  (again  in  stream  function  -  vorticity  form)  to  the 
melt  flow  field  problem.  As  in  Ref.  2,  the  melt  free  surface  was  assumed  to  be 
flat;  however,  since  the  generated  flow  field  was  the  item  of  primary  interest 
in  this  study,  the  interface  surface  was  also  assumed  flat.  The  equations  were 
solved  by  using  an  explicit  time  marching  procedure  for  the  vorticity  equation 
and  a  successive  over-relaxation  procedure  (SOR)  for  the  stream  function 
equation  at  each  time  step.  Steady  state  solutions  were  obtained  to  show  the 
effect  of  bouyancy  on  the  flow  field  for  several  demonstration  test  cases  at 
physically  unrealistic  low  Reynolds  numbers  (high  viscosities  or  small  crucibles). 

The  Langlois-Shir  effort  of  Ref.  6  was  extended  by  Langlois  in  Ref.  7  to 
cases  having  realistic  Reynolds  numbers.  Several  cases  were  considered,  and 
the  cases  run  included  bouyancy  effects,  isothermal  flows,  stationary  crucibles, 
rotating  crucibles  and  effect  of  viscosity  variation.  The  calculations  were 
run  using  the  numerical  method  described  in  Ref.  6  and  for  the  cases  considered 
this  method  required  a  very  large  amount  of  computer  run  time.  According  to 
Langlois,  each  case  of  Ref.  7  required  twenty-five  (25)  minutes  of  IBM  360/195 
computation  time  per  simulated  physical  second.  Since  the  isothermal  case  with 
rotating  crucible  was  run  for  more  than  thirty-five  (35)  simulated  seconds,  this 
case  required  over  fourteen  (14)  hours  of  run  time.  For  the  case  of  a  stationary 
crucible,  a  steady  state  was  reached  in  an  uneventful  manner.  However,  when  the 
crucible  was  rotated  in  addition  to  the  crystal,  the  secondary  flow  switched 
between  two  distinct  patterns.  In  general,  the  results  obtained  seemed  physically 
reasonable.  More  recent  studies  by  Langlois  appear  in  Refs.  8  and  9. 

The  approaches  of  Refs.  2,  6-9  demonstrate  the  ability  of  computational 
methods  to  impact  upon  the  melt  velocity  and  thermal  field  predictions. 

However,  these  efforts  were  confined  to  axisymmetric  flow  and,  therefore,  could 
not  assess  possible  three-dimensional  effects  on  the  crystal  growth  process. 
Although  the  physical  problem  is  nominally  axisymmetric,  possible  sources  of 
three-dimensional  effects  are  temperature  hot  spots  on  the  crucible  wall, 
impurities  entering  the  melt  from  the  crucible,  non-axisymmetric  crucible 
or  crystal  geometries,  misalignment  of  the  crystal  and  crucible  axes, 
misalignment  of  the  axes  from  the  vertical,  etc.  Even  when  no  such  three- 
dimensional  "forcing  functions"  are  present,  it  is  possible  that  three- 
dimensional  flow  patterns  may  occur  due  to  hydrodynamic  instabilities 
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(e.g.  see  Ref.  11).  Hence,  the  three-dimensional  aspects  of  the  problem  may  have 
significant  impact  on  the  crystal  growth  process. 

Although  the  stream  function  -  vorticity  formulation  of  Refs.  2,  6-9 
could  be  extended  to  three-dimensional  flow  in  principle,  in  practice  the  three- 
dimensional  stream  function  vorticity  approach  is  difficult,  cumbersome  and 
in  general  has  not  been  used  in  solving  three-dimensional  hydrodynamic 
problems.  A  much  more  straight-forward  and  efficient  approach  would  be 
based  upon  a  solution  of  the  Navier-Stokes  equations  with  primitive  variables, 
velocity  and  pressure  (or  density),  used  as  dependent  variables.  Such  an 
approach  is  used  here  to  demonstrate  a  model  three-dimensional  melt  flow  field 
calcualtion. 

In  addition  to  its  advantage  in  a  straight-forward  application  to  three- 
dimensional  melt  flow  field  problems,  the  velocity-pressure  approach  also 
allows  a  more  direct  analysis  to  the  free  surface  condition.  If  the  free 
surface  is  to  be  properly  modeled,  it  is  necessary  to  apply  a  pressure 
boundary  condition  there,  and  as  pointed  out  by  Langlois  (Ref.  6),  application 
of  this  pressure  boundary  condition  in  the  stream  function-vorticity 
framework  may  be  difficult.  When  primitive  variables  are  used,  application 
of  the  free  surface  boundary  condition  is  expected  to  be  considerably  easier. 
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4.  ANALYSIS 


The  crystal  growth  melt  problem  focuses  upon  thermal  and  velocity 
fields  of  a  nearly  incompressible  fluid  driven  by  inertial,  viscous  and 
bouyancy  forces.  The  bouyancy  forces  result  from  fluid  temperature  variation 
causing  a  slight  density  variation  which  leads  to  a  gravitational  body  force. 

The  present  approach  uses  the  Boussinesq  type  approximation  in  which  the 
buoyancy  force  is  modelled  as  a  body  force  proportional  to  the  difference 
between  the  local  temperature  and  a  reference  temperature.  This  approach 
has  been  used  in  Refs.  2,  6-9  and  its  validity  requires  changes  in  temperature 
to  be  small  compared  to  the  mean  temperature.  Therefore,  the  relevant  equations 
to  be  solved  are  the  incompressible  Navier-Stokes  equations  with  bouyancy 
represented  by  a  body  force  dependent  upon  the  local  temperature. 

The  present  effort  has  several  aims  including  the  demonstration  of  a 
three-dimensional  calculation  capability  and  development  of  an  efficient 
computer  cede.  A  numerical  solution  procedure  which  has  proven  to  be  efficient 
in  a  wide  range  of  fluid  mechanics  problems  is  the  linearized  block  implicit 
(LBI)  procedure  of  Briley  and  McDonald  (Refs.  12  and  13) .  Since  the  procedure 
has  been  described  in  considerable  detail  elsewhere,  it  shall  only  be  briefly 
outlined  here.  The  reader  interested  in  greater  detail  should  refer  to 
[12,13].  The  method  can  be  briefly  outlined  as  follows:  the  governing  equations 
are  replaced  by  an  implicit  time  difference  approximation,  optionally  a  backward 
difference  or  Crank-Nicolson  scheme.  Terms  involving  nonlinearities  at  the 
implicit  time  level  are  linearized  by  Taylor  expansion  in  time  about  the 
solution  at  the  known  time  level,  and  spatial  difference  approximations  are 
introduced.  The  result  is  a  system  of  multidimensional  coupled  (but  linear) 
difference  equations  for  the  dependent  variables  at  the  unknown  or  implicit 
time  level.  To  solve  these  difference  equations,  the  Douglas-Gunn  [14]  procedure 
for  generating  alternating-direction  implicit  (ADI)  schemes  as  perturbations  of 
fundamental  implicit  difference  schemes  is  introduced.  This  technique  leads  to 
systems  of  coupled  linear  difference  equations  having  narrow  block-banded 
matrix  structures  which  can  be  solved  efficiently  by  standard  block-elimination 
methods. 
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The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for  the  integration  of  initial-value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the  solution 
of  coupled  nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree 
of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no  iteration  is  required 
to  compute  the  solution  for  a  single  time  step,  and  since  only  moderate  effort 
is  required  for  solution  of  the  implicit  difference  equations,  the  method  is 
computationally  efficient;  this  efficiency  is  retained  for  multidimensional 
problems  by  using  ADI  techniques.  The  method  is  also  economical  in  terms  of 
computer  storage,  in  its  present  form  requiring  only  two  time-levels  of 
storage  for  each  dependent  variable.  Furthermore,  the  ADI  technique  reduces 
multidimensional  problems  to  sequences  of  calculations  which  are  one  dimensional 
in  the  sense  that  easily-solved  narrow  block-banded  matrices  associated  with 
one-dimensional  rows  of  grid  points  are  produced. 

One  immediate  problem  which  arises  in  the  application  of  the  LBI  procedure 
to  incompressible  problems  is  the  lack  of  a  time  derivative  of  pressure  appearing 
in  the  governing  equations.  As  shown  by  Briley  and  McDonald,  the  LBI  approach 
requires  the  appearance  of  time  derivatives  of  each  dependent  variable.  Therefore, 
if  this  procedure  is  to  be  used,  a  device  of  some  sort  must  be  created  to 
introduce  the  missing  time  derivative  into  the  equations.  One  method  of  treating 
this  problem  is  the  so-called  artificial  compressibility  approach  in  which  a 
time-derivative  of  the  required  variable  (e.g.  pressure)  is  introduced  into  the 
set  of  equations,  and  the  equations  are  marched  to  steady  state.  When  steady 
state  is  reached,  the  time  derivative  is  zero  and  the  solution  satisfies  the 
required  steady  state  equations.  Such  an  approach  has  been  used  by  several 
investigators  e.g.  Chorin  (Ref.  15). 

The  present  approach  utilizes  a  similar  4 nilosophy  in  a  somewhat 
different  context  and  is  based  upon  the  analysis  of  Briley,  McDonald  and  Shamroth 
which  shows  that  for  constant  total  temperature  flow,  the  compressible  flow 
equations  reduce  to  the  incompressible  equations  as  the  Mach  number  approaches 
zero.  This  analysis  has  been  derived  in  an  internal  SRA  report  which  is 
included  as  an  Appendix  of  the  present  report.  Following  this  analysis,  the 
present  approach  approximates  the  governing  equations  by  a  set  of  low  Mach 
number  compressible  flow  equations  and  a  gas  law  consistent  with  a  constant 
total  temperature  flow.  It  should  be  noted  that  the  constant  total  temperature 
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assumption  is  only  incorporated  into  the  assumed  gas  law  thus  making  the 
equations  of  continuity  and  momentum  approximate  the  incompressible  flow 
equations.  The  momentum  equation  contains  a  body  force  type  term  representing 
buoyancy  which  depends  upon  local  physical  temperature  and  the  local  temperature 
is  determined  from  an  energy  equation. 

With  this  introductory  paragraph  in  mind,  the  equations  used  are  the 
momentum  equation. 


,[£+(v-v>v] 


-V-P  4  F  -  V  X  |M(y  x  V)]  +  V|(*m)v-V)  (1) 


the  continuity  equation. 


ff+V(pV)  =  0 

and  the  energy  equation 


sQ 

at 


+  *  +  •v*?' 


(2) 


(3) 


where  t  represents  time,  p  represents  density,  P  represents  pressure,  V 
represents  velocity,  F  represents  body  force,  p  represents  viscosity,  e  repre¬ 
sents  internal  energy,  Q  represents  internal  heat  generation,  k  represents 
thermal  conductivity  and  4>  is  the  viscous  dissipation  function  (which  is  a 
function  of  velocity  derivatives).  D/Dt  represents  the  material  or  Stokes 
derivative.  These  equations  must  be  supplemented  by  an  equation  of  state  which 
in  general  relates  p,  p  and  T;  and  expressions  which  relate  viscosity  and  thermal 
conductivity  to  temperature. 

In  regard  to  the  equation  of  state,  the  analysis  (see  Appendix)  shows 
that  a  constant  total  temperature  fluid  approaches  an  incompressible  behavior 
as  the  Mach  number  approaches  zero.  Therefore,  the  equation  of  state  is 
taken  as 

p'-fR(T~^)  <4> 


where  R  is  the  gas  constant  and  is  specific  heat.  If  t  were  the  local 
total  temperature  in  the  fluid,  Eq.  (4)  would  represent  the  perfect  gas  law 
for  an  ideal  fluid.  However,  in  the  present  case  t  is  taken  as  a  ficticous 
constant  temperature  of  high  enough  value  so  that  the  Mach  number,  based 
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upon  V  and  t  is  small;  i.e., 

Mf  =  1VI  / yVR  (t-  V2/2Cp)  «  1.0 

In  the  present  calculations  this  Mach  number  was  kept  below  0.03  and  consequently 
the  equations  approach  their  incompressible  form.  The  final  item  to  be 
considered  is  the  bouyancy  force.  In  the  present  calculations,  the  bouyancy 
force  is  represented  by  the  Boussinesq  relation 

7= +^>g'/3(T-Tr)  (5) 

where  6  is  the  volumetric  coefficient  of  expansion,  g  is  the  acceleration  due 
to  gravity,  T  is  the  local  temperature  and  is  a  reference  temperature. 

The  boundary  conditions  used  in  the  present  analysis  consist  of  no-slip 
and  one-sided  momentum  equations  on  all  solid  surfaces.  On  the  free  surface, 
the  velocity  normal  to  the  surface  is  aero,  the  normal  derivatives  of  the  other 
velocity  components  are  zero  and  a  one-sided  momentum  equation  is  solved.  The 
temperature  is  specified  on  both  the  crystal  and  crucible  surfaces  and  a 
second  normal  derivative  of  temperature  is  set  to  zero  on  the  free  surface. 

The  equations  are  solved  as  a  coupled  system  with  convergence  acceleration 
techniques  of  the  type  described  by  McDonald  and  Briley  (Ref.  16). 
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5.  RESULTS 


Several  calculations  were  run  under  the  present  effort.  These  included 
a  code  verification  test  case,  axisynnnetric  cases  with  and  without  bouyancy 
effects  and  a  three-dimensional  demonstration  case.  These  are  now  discussed 
in  detail. 

The  first  case  was  run  to  verify  the  code  against  a  known  analytic 

solution.  The  case  considered  was  flow  within  a  rotating  crucible  with  no 

crystal  in  the  free  surface.  For  this  case  the  analytic  solution  for  the 

azimuthal  velocity  component  is 

V  (r)  =  V  - 
0  w  r 

w 

where  V  is  the  velocity  at  the  wall  and  r  is  the  wall  radius.  The  calculation 
w  w 

was  run  by  assuming  an  initial  flow  and  marching  the  equations  in  time  until 
convergence  was  obtained.  Convergence  was  judged  by  monitoring  the  maximum 
flow  field  residual  as  well  as  the  maximum  change  in  dependent  variables  over 
a  time  step.  A  comparison  between  the  predicted  V.  profile  and  that  obtained 
from  the  analytic  solution  for  an  incompressible  fluid  is  given  in  Table  I. 

As  can  be  seen,  the  numerical  solution  is  in  excellent  agreement  with  analysis. 
The  calculation  was  repeated  for  flow  between  co-rotating  cylinders  and  the 
comparison  between  the  numerical  predictions  and  the  known  analytic  solution 
again  was  very  good.  These  comparisons  serve  to  verify  the  code  for  these 


simple  test  cases. 


The  second  case  considered  was  axisymmetric  flow  in  a 

stationary 

crucible  with  the  crystal  rotating. 

the  crucible 

If  the  crystal  radius 

ve 

is  denoted  by  R 

r 

Predicted 

Analytical 

1.0 

1.000 

1.0 

.95 

.9499 

.95 

.90 

.8999 

.90 

.85 

.8499 

.85 

.80 

.7999 

.80 

.75 

.7499 

.75 

.70 

.6989 

.70 

.65 

.6499 

.65 

.60 

.5988 

.60 

.55 

.5498 

.55 

.50 

.4998 

.50 

.45 

.4498 

.45 

.40 

.3998 

.40 

.35 

.3498 

.35 

.30 

.2998 

.30 

.25 

.2499 

.25 

.20 

.1999 

.20 

.15 

.1499 

.15 

.10 

.0999 

.10 

.05 

.0499 

.05 

.0 

.0 

.0 

Table  1  -  Comparison  of  predicted  and  analytical  velocity  profiles  for  rotating  crucibles 
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The  Reynolds  number  Re  represents  the  ratio  of  inertial  to  viscous  forces  and 

this  value  represents  a  typical  crystal  growth  process  in  which  R.  =  4.2cm, 

1  2 

R^  =  6.4cm,  V1  =  9.7cm/sec,  =  0  lOcm/sec ,  h  =  3.95cm  and  v  =  0.01cm  /sec. 
Bouyancy  effects  were  set  to  zero  for  this  case.  With  these  parameters  the 
case  corresponds  to  the  stationary  crucible  case  considered  by  Langlois  in  Ref.  6. 

The  calculation  run  on  a  grid  containing  thirty-one  (31)  radial  points  and 
twenty-seven  (27)  axial  points  converged  to  a  steady  solution  within  eighty  (80) 
time  steps.  The  present  deck  is  a  very  general  code  designed  for  maximum 
flexibility  in  which  equations,  dependent  variables,  grid  specification, 
boundary  conditions,  etc.  can  be  easily  changed  and,  therefore,  the  code 
contains  a  great  deal  of  computer  overhead.  Even  with  this  large  amount  of 
overhead,  the  code  required  less  than  five  (5)  seconds  of  CDC  7600  time  per 
time  step.  Recent  changes  designed  to  decrease  run  time  (accompanied  by  a 
decrease  in  code  flexibility)  have  decreased  run  time  by  a  factor  of  2.  If 
the  calculation  were  repeated  with  this  new  code,  the  total  run  time  would  be 
approximately  3  minutes.  It  is  anticipated  that  further  streamlining  of  the 
code  will  reduce  run  time  by  an  additional  factor  of  3. 

Predictions  of  the  velocity  and  temperature  fields  are  shown  in  Figs.  2-4. 
Both  the  Langlois  calculation  and  the  present  calcualtion  reached  a  steady 
state.  The  azimuthal  velocity.  Fig.  2,  shows  good  qualitative  agreement 
with  the  Langlois  solution.  Both  calculations  show  the  thin  boundary  layer 
on  the  crystal  face  as  well  as  the  local  maximum  near  the  stationary  crucible 
wall.  The  predicted  secondary  flow  field  (Figure  3)  shows  the  velocity 
components  and  is  also  in  good  qualitative  agreement  with  the  results 

of  Langlois.  Both  calculations  show  the  secondary  flow  to  consist  of  a  single 
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vortex  pattern  which  is  centered  in  the  upper  right  corner  of  the  flow  field. 
Finally,  temperature  contours  are  shown  in  Fig.  4  where  the  ratio  of  the 
crucible  wall  temperature  to  the  crystal  face  temperature  is  1.052. 

The  next  case  considered  corresponds  to  Langlois  rotating  crucible  case 
without  buoyancy  forces.  The  parameters  were  now  changed  so  that 

V2/Vl  =  -1*02 

All  other  parameters  correspond  to  that  of  the  first  case  and  again  bouyancy  is 
neglected.  When  bouyancy  is  neglected  the  velocity  field  is  unaffected  by 
the  thermal  field  and  hence  the  resulting  velocity  field  represents 
an  isothermal  calculation.  The  calculation  was  again  run  on  a  31  x  27 
computational  grid,  however,  in  this  case  complete  convergence  was  not  obtained. 

The  major  region  of  oscillation  occurred  near  the  free  surface  at  radii  just 
beyond  the  crystal  boundary,  (R^) .  The  velocity  in  this  region  was  not  well 
resolved  and  it  is  likely  that  this  lack  of  resolution  may  have  prevented 
convergence.  However,  the  calculation  did  not  show  the  unsteady  flow  pattern 
observed  in  Langlois  case  where  the  secondary  flow  shifted  between  two 
distinct  patterns. 

The  azimuthal  velocity  contours  for  this  case  are  shown  in  Fig.  5.  The 
contours  show  the  Taylor  column  under  the  crystal  face  and  are  in  good 
qualitative  agreement  with  the  results  of  Langlois.  The  four  different 
regions  within  the  Taylor  column  discussed  by  Langlois  are  clearly  evident. 

Directly  under  the  crystal  fluid  turns  with  the  crystals  rather  than  with 
the  crucible.  This  is  followed  by  a  transition  region  in  which  the  flow 
azimuthal  velocity  changes  from  the  crystal  to  the  crucible  direction.  Next 
comes  a  region  where  the  velocity  changes  little  with  depth  and  finally  the 
boundary  region  at  the  crucible  bottom. 

The  secondary  flow  vector  plot  in  Fig.  6  shows  four  distinct  vortices: 
a  strong  clockwise  vortex  just  under  the  crystal  face,  a  large  but  weaker 
counterclockwise  vortex  encompassing  most  of  the  region  between  the  first 
vortex  and  the  crucible  bottom  and  two  weak  vortices  in  the  outer  radial  region. 

These  latter  two  vortices  are  very  much  elongated  in  the  z-direction.  Langlois* 
solution  showed  a  weak  vortex  in  the  outer  region  which  shifted  between  two  patterns. 
The  maximum  values  of  the  predicted-secondary  flow  velocities  are  approximately 
0.10  within  the  Taylor  column  and  0.02  outside  of  the  Taylor  column. 
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These  values  are  in  good  agreement  with  those  given  by  Langlois.  Temperature 
profiles  are  presented  in  Fig.  7. 

The  final  axisymmetric  case  considered  adds  bouyancy  effects  to  the 
calculation.  The  flow  parameters  used  were  those  of  the  proceeding  case  with 
the  following  additions 

Pr  =  0.08 

—  (T  -  T  )  r  1  0 

^crucible  crystal' 

The  results  are  presented  in  Figs.  8-10.  The  azimuthal  velocity  field  shown 
in  Fig.  8  again  is  in  qualitative  agreement  with  Langlois'  calculation.  In 
particular  the  change  due  to  bouyancy  effects  (see  Figs.  5  and  8)  correspond 
to  the  changes  shown  by  Langlois.  The  most  dramatic  change  due  to  bouyancy 
occurs  in  the  secondary  flow  profile  where  as  shown  in  Fig.  9,  the  strength 
of  the  secondary  flow  pattern  is  increased  considerably.  The  temperature 
field  in  this  case  shown  in  Fig.  10.  These  also  show  good  qualitative  agreement 
with  Langlois. 
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6.  THREE  DIMENSIONAL  CALCULATION 

The  final  calculation  focused  upon  three  dimensional  effects.  Before 
discussing  this  calculation  several  points  should  be  considered.  First  of  all, 
three-dimensional  calculations  make  major  demands  on  computer  resources  and, 
therefore,  for  this  demonstration  calculation  a  flow  field  was  chosen  to  allow 
the  calculation  to  be  made  with  a  moderate  number  of  grid  points  (an  11  x  17  x  16 
grid  was  used) .  This  relatively  small  number  of  three-dimensional  grid  points 
required  the  Reynolds  number  to  be  smaller  than  would  normally  be  expected  in  a 
physical  problem.  Therefore,  the  case  must  be  viewed  as  a  demonstration  case. 

The  second  point  concerns  the  source  of  the  three  dimensionality.  It 
was  originally  anticipated  that  the  three  dimensionality  would  result  from  the 
angular  misalignment  of  the  crystal  and  crucible.  The  present  effort  assumes 
a  flat  free  surface  and,  therefore,  the  effect  of  angular  misalignment  would  be 
confined  to  changing  the  region  on  the  upper  surface  bounded  by  the  crystal 
from  a  circle  to  an  ellipse.  For  a  5-degree  misalignment,  the  resulting 
ellipse  would  have  a  ratio  of  major  to  minor  axis  of  1.004.  In  the  absence 
of  meniscus  effects,  the  three  dimensionality  arising  from  misalignment  would 
be  negligible  and,  therefore,  a  second  three  dimensional  source  arising  from  a 
hot  spot  on  the  crucible  wall  was  added  to  the  calculation. 

The  calculation  was  run  using  16  radial  grid  points,  11  axial  grid  points 
and  17  azimuthal  grid  points.  The  crystal  was  assumed  to  be  rotating  and  the 
crucible  stationary. 

The  geometrical  parameters  were  taken  as 

R2/R1  =1.5 
h/R1  =0.5 

The  flow  parameters  were  taken  as 

Re  =  V^/v  =  100 

Pr  =  .05 

£ 

Ff  ^crubile  -  ^crystal)  =  .02 

where  is  the  velocity  at  the  crystal  radius,  R^,  v  is  viscosity,  Pr  is 
Prandtl  number,  B  is  the  volumetric  coefficient  of  expansion  and  Fr  is  the 
Froude  number.  The  ratio  of  the  crucible  wall  temperature  to  the  crystal 
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temperature  was  taken  to  be  1.052  and  the  hot  spot  temperature  ratio  was 
taken  as  1.157.  In  order  to  further  limit  the  computational  domain  two 
identical  hot  spots  were  assumed  present.  These  were  centered  at  z/h  =  0.4, 

8  =  n/2  and  z/h  =  0.4,  8=311/2.  This  specification  allowed  the  calculation 
to  be  confined  to  the  domain  0  <  8  <  n. 

The  calculation  was  first  converged  as  an  axisymmetric  calculation  with 
no  axis  angular  misalignment  and  constant  crucible  wall  temperature.  The 
azimuthal  velocity  contours,  secondary  flow  vector  plots  and  temperature  contours 
are  shown  in  Figs.  11-13.  The  results  are  as  expected  with  the  low  Reynolds  number 
allowing  significant  azimuthal  velocity  to  occur  through  most  of  the  flow  field 
(compare  this  to  Fig.  2  where  Re  =  4000).  A  wall  hot  spot  was  then  added  to 
the  problem  and  a  three-dimensional  calculation  initiated  from  the  converged 
two-dimensional  calculations.  Although  both  the  velocity  and  temperature 
fields  showed  a  three-dimensional  variation,  the  velocity  field  variation  in 
the  azimuthal  direction  was  very  small.  Both  the  azimuthal  velocity  contours 
and  secondary  flow  pattern  corresponded  very  closely  to  the  axisymmetric  case. 

The  temperature  field  showed  significant  variation  and  the  temperature  field 
in  various  azimuthal  planes  are  shown  in  Figs.  14-16.  The  calculation  then  was 
continued  with  a  5-degree  angular  misalignment  between  the  crystal  and  crucible 
axes;  as  previously  discussed,  the  addition  of  this  misalignment  factor  was 
not  large  enough  to  change  the  predicted  flow  field. 
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7.  CONCLUSIONS 

The  present  effort  has  applied  a  Navier-Stokes  analysis  using  velocity  and 
density  as  dependent  variables  to  the  crystal  growth  melt  problem..  The  approach 
models  the  equations  governing  the  velocity  field  as  low  Mach  number  compressible 
flow  equations  with  a  body  force  proportional  to  the  local  temperature  and  a  gas 
law  which  approximates  incompressible  flow.  The  analysis  has  been  applied  to 
flow  within  a  rotating  cylinder  and  flow  between  two  concentric  cylinders,  and 
in  both  cases  the  predicted  azimuthal  velocity  field  showed  excellent  agreement 
with  the  analytic  solution  to  the  governing  equations.  Axisymmetric  melt 
problems  which  included  the  effects  of  bouyancy  and  crucible  rotation  were 
run  and  the  predicted  velocity  and  thermal  fields  showed  good  qualitative 
agreement  with  the  vorticity-stream  function  solutions  of  Langlois.  The 
approach  was  then  applied  to  a  model  three-dimensional  problem  in  which  the 
three  dimensionality  resulted  from  a  hot  spot  on  the  crucible  wall.  Although 
the  flow  parameters  (particularly  the  low  Reynolds  number)  required  this  to  be 
viewed  as  a  model  calculation,  the  calculation  definitely  showed  the  capability 
of  the  procedure  to  calculate  three-dimensional  flows.  An  assessment  of  the 
effects  of  three-dimensional  sources  upon  the  crystal  growth  process  will 
require  further  code  evaluation,  code  'speed-up'  and  extension  of  the  code 
to  include  items  such  as  the  free  surface  boundary. 
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Fig.  3  -  Secondary  flow,  stationary  crucible  cas 
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Pig.  8  -  Azimuthal  velocity  contours,  rotating  crucible  with  bouyancy. 
Contour  levels  at  -.8,  -.5,  -.2,  0,  .1,  .2,  .5,  .75,  .9. 


Pig.  9  -  Secondary  flow,  rotating  crucible  with  buoyancy. 


Fig.  14  -  Temperature  contours  at  0  ■  36®,  three-dimensional  case  (hot  spot). 

Contour  levels  at  1.00.  1.01.  1.02.  1.03.  1.04.  1.05.  1.06.  1.07.  1 


Temperature  contours  at  0  ■  84* ,  three-dimensional  case  (hot  spot). 
Contour  levels  at  1.00,  1.01,  1.02,  1.03,  1.04,  1.05,  1.06,  1.07,  1 
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ABSTRACT 

It  is  well  known  that  flow  of  a  perfect  gas  at  low  Mach  number  closely 
approximates  an  incompressible  flow,  provided  there  is  no  heat  addition.  The 
usual  formulation  of  the  compressible  Euler  equations  encounters  singular  behavior 
as  the  Mach  number  approaches  zero.  The  Euler  equations  for  a  perfect  gas  are 
written  here  in  variables  which  are  well  behaved  for  adiabatic  flow  at  low  Mach 
number  and  in  a  form  which  reduces  to  the  incompressible  (constant  density) 

Euler  equations  as  the  Mach  number  approaches  zero.  The  present  formulation 
clarifies  how  the  reduction  to  incompressible  flow  occurs  and  is  useful  In 
constructing  numerical  algorithms  for  low  Mach  number  and  other  flows. 
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BACKGROUND 


Incompressible  Euler  Equations 

The  incompressible  Euler  Equations  can  be  expressed  in  nondimensional 


vector  form  as 


7  •  U  =  0  (1) 

||  +  <0  •  V)  0  +  V  =  0  (2) 

where  U  is  the  velocity  vector,  c^  is  a  pressure  coefficient,  and  t  is  time. 

The  pressure  coefficient  is  defined  by 

cp  =  2  (p  -  Px)  (3) 

where  p  is  pressure  and  P^  is  an  arbitrary  (constant)  pressure  basis.  Here 
(and  subsequently),  all  variables  are  nondimensional,  having  been  normalized  by 

reference  quantities  denoted  by  a  subscript  'r'.  The  reference  pressure  p  has 

2  r 
been  taken  as  p^U^  «  where  p  denotes  density. 

The  continuity  equation  (1)  and  the  three  components  of  momentum  (2) 

comprise  four  equations  governing  c^  (or  equivalently  p)  and  the  three 

components  of  U.  For  incompressible  flow,  c^  and  U  are  independent  of  the 

temperature  T,  which  in  turn  is  governed  by  an  energy  equation  that  does  depend 

on  U. 

Compressible  Euler  Equations 

The  compressible  continuity  and  momentum  equations  can  be  written  in 
the  following  form: 


9(£n  p) 


+  V  •  U  +  U  •  V  (£n  p)  =  0 


1“  +  U  +  (u  •  V) 


u  +  V  +  £  V  Un  p)  =  0  (5) 


The  equation  of  state  for  a  perfect  gas  can  be  expressed  as 

P  =  pT/yM2  (6) 

where  y  is  the  specific  heat  ratio,  and  the  reference  Mach  number  M  is  given  by 

2 

M  =  Ur/c.  Here,  c  is  the  reference  speed  of  sound  defined  by  c  =  yRT^.,  and  R 
is  the  gas  constant.  The  'r'  svibscript  on  M  is  omitted  for  convenience,  but  no 
confusion  should  arise  since  only  the  reference  Mach  number  will  be  referred  to 
subsequently.  The  system  is  closed  by  the  energy  equation,  which  may  be 
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written  as 


(7) 


+  0  •  VT  =  (y-l)T  V  •  U 
3 1 

The  usual  method  of  reducing  the  compressible  equations  to  the  incompres¬ 
sible  equations  is  to  assume  a  constant  density  and  note  that  an  equation  of 
state  is  no  longer  needed.  As  a  limiting  process,  however,  the  equation  of 
state  (6)  indicates  that  p  becomes  infinite  as  M  -*■  0  (unless  T  0)  .  Equations 
(4  -  7)  are  reformulated  here  as  a  system  of  equations  which  remains  well 
behaved  as  M  ->  0. 


LOW  MACH  NUMBER  FORMULATION 


To  relate  the  compressible  and  incompressible  equations,  a  new 
dependent  variable  c^  relating  p  and  p  is  introduced  and  defined  by 


c  = 

P 


yM' 


2  j  W2*  '  V 


(8) 


where  P ^  is  an  arbitrary  (constant)  pressure  basis.  Using  Eq.  (8)  to  eliminate 
p  in  the  momentum  equation  (5)  leads  to 

A  A 

_  c  c 

U  +  u  +  (u  •  V)U  +  V  +-^  V(ta  p)  »  0  (9) 

The  definition  of  (nondimensional)  stagnation  enthalpy  E  is 

E  =  T  +  (y-1)  M2  q2/2  (10) 

2  -  - 

where  q  =  U  •  U,  and  E  has  been  taken  as  C  T  ,  where  C  is  the  specific 
^  r  p  r  p 

heat  at  constant  pressure.  Combining  the  equation  of  state  (6)  and  the  defini¬ 
tion  of  stagnation  enthalpy  (10)  gives 

p  =  p  [E/yM2  -  (y-l)  q2/2y]  (11) 

Combining  Eqs.  (8)  and  (11)  to  eliminate  p  leads  to 

P  -  p2  it  -  ^  p2]-1  a?) 

For  a  wide  range  of  flow  conditions,  including  adiabatic  flow  at  low  Mach 
number,  the  stagnation  enthalpy  E  may  be  assumed  to  be  a  known  constant.  In  this 
instance,  Eq.  (12)  relates  p,  c  ,  and  U  and  serves  to  decouple  the  continuity 
and  momentum  equations  (4)  and  (9)  from  the  energy  equation  (7),  which  can  be 
omitted  from  consideration.  The  system  of  governing  equations  then  consists  of 
Eqs.  (4),  (9),  and  (12),  with  dependent  variables  p,  cp,  and  U. 
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Note  from  Eq.  (12)  that  p  P2/E  (a  constant)  as  M  +  0.  As  a  consequence, 
the  terms  containing  V(£n  p)  in  Eqs.  (4)  and  (9)  vanish  as  M  -*■  0,  and  these 
equations  thus  reduce  to  the  incompressible  form  of  Eqs.  (1)  and  (2).  The  limit 
M  -*■  0  is  taken  by  letting  R  °°,  which  allows  U^.,  pr,  T  ,  y  =  0(1). 

In  this  formulation, 

c  c  ,  p  1  as  M  +  0  (13) 

P  P 

assuming  appropriate  choices  are  made  for  P  ,  P„ ,  and  E.  To  guide  these  choices, 

it  is  noted  that  selecting  =  1,  E  =  1  +  (y-l)M  /2  implies  that  c^  =  0  and 

p  =  T  =  1  where  q  =  1  (freestream  basis  condition  for  P^).  The  choice  =  E  =  1 

implies  that  c^=  0  and  p  =  T  =  1  where  q  =  0  (stagnation-point  basis  for  P^). 

With  either  choice  for  P_  and  E,  c  c  follows  as  M  0  because  c  and  c  satisfy 

2  p  p  P  P 

the  same  governing  equations  and  boundary  conditions  and  are  equal  at  the  basis 
point.  Also,  since  p  -*•  1  as  M  -*  0,  £n  p  itself  approaches  zero  as  M  0. 

It  is  important  to  distinguish  the  roles  played  by  c^and  by  the  assumption 
of  constant  E.  The  change  of  variables  from  p  to  cp  is  responsible  for  eliminat¬ 
ing  the  singular  behavior  which  would  otherwise  occur  in  Eq.  (11)  whether  or  not 
E  is  constant.  The  assumption  of  constant  stagnation  enthalpy  is  not  a  requirement 
for  use  of  the  low  Mach  number  formulation  of  the  equations,  and  is  necessary  only 
if  the  compressible  equations  are  to  imply  constant  density  as  M  -*■  0  (cf.  Eq.  12). 
It  is  thus  clear  that  the  physical  motivation  for  requiring  constant  stagnation 
enthalpy  at  low  Mach  number  is  to  produce  a  flow  field  with  nearly  constant 
density. 

It  is  also  worth  noting  that  since  V  •  D  +  0  as  M  +  0,  the  last  term  in 
the  energy  equation  (7)  vanishes,  and  this  equation  thus  reduces  to  its  incompres¬ 
sible  form 

~  +  U  •  VT  =  0  (14) 

O  t 

Since  Eqs.  (10)  and  (12)  imply  that  T  -*  E  and  p  P^/E  as  M  0,  it  is  evident 
that  taking  E  as  constant  as  M  +  0  does  not  permit  a  variable  T,  and  if  E  varies, 
then  p  varies  as  T-^.  Nevertheless,  the  present  formulation  can  be  used  to  treat 
incompressible  (constant  density)  flows  with  nonconstant  temperature  by  taking  E 
as  constant  and  omitting  Eq.  (10).  This  can  be  justified  by  noting  that  for  small 
Mach  number  (say  M  <  0.1)  and  constant  E,  the  present  system  for  a  perfect  gas 
closely  approximates  a  constant  density  incompressible  flow  and  is  thus  a  good 
model  (for  example)  for  flow  of  a  liquid  whose  (constant)  density  is  assumed  to  be 
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independent  of  temperature.  In  this  instance,  the  quantity  E  in  Eq.  (12)  is 
reinterpreted  as  a  nonphysical  modeling  parameter  which  serves  to  decouple  the 
continuity  and  momentum  equations  from  the  energy  equation,  as  required  for 
the  liquid.  The  incompressible  energy  equation  (14)  can  then  be  solved  separately 
for  the  physical  temperature  variation,  and  the  incompressible  constant-density 
liquid  is  thus  modeled  as  a  very  slightly  compressible  fluid  with  constant 
'stagnation  enthalpy'.  The  Boussinesq  assumption  for  natural  convection  can  be 
introduced  within  this  same  context  by  adding  a  body  force  term  (which  depends 
only  on  temperature)  to  the  momentum  equation. 

As  a  separate  issue,  it  is  worth  noting  that  the  assumption  of  constant 
stagnation  enthalpy  is  useful  over  a  wide  range  of  flow  conditions  and  is  much 
more  general  than  is  an  assumption  of  isentropic  flow.  The  constant  stagnation 
enthalpy  assumption  is  reasonable  for  transonic  and  supersonic  flow  with  shocks 
and  also  for  adiabatic  viscous  flow  with  unit  Prandtl  number  and  no  viscous 
heating.  The  authors  have  made  extensive  use  of  this  assumption  in  treating 
viscous  flow  at  supersonic,  low  subsonic,  and  transonic  Mach  numbers  (eg.  Refs.  1-3). 
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